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Abstract 


Disturbance events strongly affect the composition, structure, and function of forest 
ecosystems; however, existing U.S. land management inventories were not designed to 
monitor disturbance. To begin addressing this gap, the North American Forest Dynamics 
(NAFD] project has examined a geographic sample of 50 Landsat satellite image time series 
to assess trends in forest disturbance across the conterminous United States for 1985- 
2005. The geographic sample design used a probability-based scheme to encompass major 
forest types and maximize geographic dispersion. For each sample location disturbance 
was identified in the Landsat series using the Vegetation Change Tracker (VCT] algorithm. 
The NAFD analysis indicates that, on average, 2.77 Mha/yr of forests were disturbed 
annually, representing 1.09%/yr of US forestland. These satellite-based national 
disturbance rates estimates tend to be lower than those derived from land management 
inventories, reflecting both methodological and definitional differences. In particular the 
VCT approach used with a biennial time step has limited sensitivity to low-intensity 
disturbances. Unlike prior satellite studies, our biennial forest disturbance rates vary by 
nearly a factor of two between high and low years. High western US disturbance rates were 
associated with active fire years and insect activity, while variability in the east is more 
strongly related to harvest rates in managed forests. We note that generating a geographic 
sample based on representing forest type and variability may be problematic since the 


spatial pattern of disturbance does not necessarily correlate with forest type. We also find 
that the prevalence of diffuse, non-stand clearing disturbance in US forests makes the 
application of a biennial geographic sample problematic. Future satellite-based studies of 
disturbance at regional and national scales should focus on wall-to-wall analyses with 
annual time step for improved accuracy. 
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1.0 Introduction 

Change is ubiquitous in forest ecosystems. Forests experience both seasonality as well 
as long-term growth cycles that can vary in duration between 50 years and 500 or more 
years (Waring & Running, 2007], These long-term changes are punctuated by mostly short- 
term disturbances from fire, insects, disease, and harvest which strongly alter the state and 
functioning of the forest (He & Mlandenoff, 1999]. Both climate change and the increasing 
global demand for wood and fiber products are likely to drive increases in forest 
disturbance rates (Kurz et al., 2008; Nepstad et al., 2008], These changes in disturbance will 
alter the water and carbon cycles of forest stands as well as impact the habitat and 
biodiversity of these ecosystems (Lindenmayer et ah, 2006; Gardner et ah, 2009]. With 
respect to the carbon cycle, forest disturbance is now recognized as a major driver of non- 
fossil-fuels-related terrestrial fluxes to the atmosphere (Running, 2008; Amiro et ah, 2010]. 

To effectively understand how forest disturbance impacts forest state and functioning, 
disturbance rates need to be quantified at the spatial grain where human management and 
natural disturbances occur; typically less than 10 ha (Miller, 1978; Cohen et ah, 2002; 
Kuemmerle et ah, 2007; Frolking et ah, 2009]. Further, disturbances need to be quantified 
at a time step relevant to the various affected processes (e.g., annually or more frequently], 
and at spatial (national to global] and temporal (e.g. historically back at least as far as the 
1990s] extents relevant to establishment of baselines meaningful to forest policy initiatives 
(Bottcher et ah 2008, Masek et ah 2008, Kennedy et ah 2012]. 
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In the United States, the lack of consistent, high-temporal-resolution estimates of forest 
disturbance remains an important gap in efforts to model and manage forest carbon at a 
national scale (USCCSP, 2007; Birdsey et al., 2009]. The US Forest Service Forest Inventory 
and Analysis (FIA] Program relies on a network of plots to inventory and monitor forested 
ecosystems at regional to national scales. Its current annual inventory system (McRoberts 
et al. 2005] closely tracks individual tree mortality and cause of disturbance through 
remeasurement of inventory plots on a regular cycle (approximately every 5 years in the 
east, ten in the west], ffowever, until remeasured data is available nationally, consistent 
forest disturbance estimates cannot be constructed. In addition, FIA is not structured to 
capture relatively rare disturbance events. Consequently, today’s reported US national 
inventory-based estimates of disturbance area have drawn upon separate databases for the 
extent of harvest (Smith et al., 2009], fire (US EPA, 2011], and insect damage (USDA Forest 
Service, 2010], In some cases these estimates may be inconsistent. For example, harvest 
area is derived from a combination of inventory re-measurement data, where available, and 
harvest activity reports from different National Forests. Insect mortality typically reflects 
the gross area affected by insects as measured by aerial sketch maps over purposively 
sampled regions of the country (Johnson and Wittwer 2008], although recent efforts have 
begun to convert these maps into true mortality estimates (Meddens et al., 2012], 

Satellite observations may provide a more consistent means for assessing disturbance. 
Previous studies estimating national and global forest disturbance patterns have used 
coarse-resolution (250m - 1km] satellite imagery (NOAA AVHRR and NASA EOS MODIS] 
(Potter et al., 2005; Mildrexler et al., 2009; Potapov et al., 2009]. Although suitable for 
detection of large-area disturbances such as fire and large-scale clearcuts, coarse-resolution 
imagery is less capable of detecting forest management activities than finer-grained 
observations from satellites such as Landsat (Skole and Tucker, 1993; Tucker and 


4 


Townshend, 2000; Bucha & Stibig, 2008; Wulder et al., 2008; Potapov et al., 2009]. Studies 
that use imagery of finer spatial grain have either mapped at too coarse a temporal grain to 
detect short-term changes in forest disturbance rate (Masek et al., 2008; Hansen et al., 

2010] or have focused only a single type of forest disturbance (e.g. fire mapping from the 
Monitoring Trends in Burn Severity project (MTBS, Eidenshink et al., 2007], The need to 
overcome limitations inherent in these previous or ongoing ground-based and satellite- 
based disturbance analysis efforts has motivated the development of the North American 
Forest Dynamics (NAFD] project (Goward et al., 2008], a core project of the North American 
Carbon Program (Wofsy & Harris, 2002]. 

The study reported here was derived from the first two phases of NAFD, which 
employed a geographic sample of Landsat observations at a relatively high temporal 
frequency (approximately annual time step] over a 20-year period to characterize the 
dynamics of recent US forest disturbance history. A sampling approach was selected due to 
the high cost of individual Landsat images (~$600] when this study was initiated (2005] 
and the fact that an annual, wall-to-wall analysis would involve over 9000 such images 
covering the 442 individual scenes over the conterminous US 1 . We selected a sample of 50 
scenes and for each we assembled a time series of cloud-minimized, seasonally consistent 
imagery. The time series of each pixel in a given image stack was analyzed using the 
Vegetation Change Tracker (VCT] algorithm to detect disturbances (Huang et al., 2010a]. 

In this study, forest disturbance was defined as any event that caused either substantial 
mortality or leaf-area reduction within a forest stand, including management activities such 
as harvest and thinning. As described more fully in the methods section, the VCT approach 
captures most rapid stand-clearing events (including clearcut harvests and fire], as well as 

1 In this paper we use the term "scene" to refer to the nominal geographic area of a 
particular Landsat Worldwide Reference System [ WRS-2 path/row). The actual Landsat 
acquisitions from specific dates for that area are termed "images" (see Strahler et al 1986). 
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many non-stand-clearing events (partial harvest, thinning, storm damage, insect damage). 
However, gradual declines in live biomass that occurred over several years (e.g. due to 
drought or disease) were mostly not captured. The approach also did not distinguish 
between disturbance (mortality followed by recovery) and permanent conversion of land 
cover. Thus our definition of disturbance corresponds most closely to "gross forest cover 
loss" (Hansen, 2010). Although knowing the causal agent of disturbance is extremely 
important for understanding specific impacts on ecosystems, this study has not attempted 
to assign an agent to each disturbed patch. Instead we focus on overall "turnover” of live 
forest area across the county. 

Our central objective was to estimate annual rates of forest disturbance across the 
conterminous US. Accomplishing this required the development and application of novel 
methods to: (i) select a probability-based sample of 50 scenes while satisfying diverse 
analytical criteria; (ii) apply an automated change detection algorithm to identify forest 
disturbance across a Landsat image time series for each scene, and (iii) assign estimates of 
sampling errors associated with this disturbance mapping. 

2.0 Methods 

2.1. Sample Design 

Our sample design followed the rules of probability-based sampling (Sarndal et al., 
1992), where each scene had a known, non-zero, positive probability of inclusion in the 
sample, and the set of possible samples was finite and known. This approach allowed for 
design-based estimation and for preferential inclusion of scenes with greater land or forest 
area (Gallego, 2005), an important consideration given the effort and cost involved in the 
analyses and the central goal of characterizing forest disturbance. This approach also 
allowed for inclusion of other important characteristics: spatial dispersal of scenes to 
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minimize autocorrelation, assurance that all major forest types were included as forest 
disturbance and recovery dynamics are largely type-specific, and the ability to leverage 
work already completed at a handful of "targeted" scenes available from related projects 
(e.g. Masek and Collatz 2006, Eidenshink et al., 2007]. 

Fundamental to our sampling design was the choice of sample unit (geographic area of a 
single sample scene] and sample frame (the population of sample scenes covering the 
conterminous US]. Because adjacent WRS-2 Landsat frames overlap, the sample frame was 
modified such that each scene was trimmed to include only the unique, non-overlapping 
area it contained. This ensured that each sample scene was unique, and simplified the use 
of population statistics. Hereafter, our reference to sample scenes implies the non- 
overlapping portion of scenes produced by this process. 

As a first step in assuring that scenes with greater forest were preferentially selected as 
part of the sample, several scenes with extremely low forest cover were culled from the 
sample frame. This was accomplished by ranking each sample by total forest area (based on 
a US forest type map developed by Ruefenacht et al, 2008] from low to high, and eliminating 
all ranked scenes for which the total forest area was below 2% of the cumulative forest area 
across all scenes in the frame. 

We divided the culled sampling frame into two strata: eastern and western. This 
division reflects that fact that at a gross scale eastern and western US forests are 
fundamentally different ecologically and in terms of how they are managed, both of which 
could greatly affect disturbance and recovery dynamics. To accomplish this division all 
Landsat scenes in WRS-2 Path 31 or greater were declared western and all scenes in Path 
30 or less were declared eastern. All further steps were executed separately for each 
stratum. 
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For each stratum, samples were drawn from randomly-ordered lists of scenes. Rather 
than creating a single list, a set of v (=100,000) scene lists were created, each with a 
different random ordering, and the final sample was chosen by randomly selecting one list 
from the set of lists that met our preferential criteria: geographic scene dispersion, 
maximizing total forest area, forest type diversity, and inclusion of targeted scenes where 
we have important experience. This provides two important benefits. First, for any given set 
of lists, every scene's probability of inclusion (jt) in the sample could be calculated directly 
as its proportional occurrence in the first n scenes across all lists in the set. Second, because 
the probability of inclusion was calculated directly, it was possible to further cull the set of 
lists from the frame to remove those whose first n scenes, taken together, did not meet our 
preferential selection as described above. 

The final sample was chosen by randomly selecting one list per stratum from the final 
culled list in each stratum. The sample size (n) for each stratum was 12 (eastern) and 11 
(western), given that the sample design and list selection was executed during the first 
phase of the project when only 23 scenes were under consideration. For each sample scene 
jt was calculated as the number of times the sample appeared among the first n scenes in 
the set of v lists, divided by v. During the second phase of the project, when we added the 
next sequential 13 (eastern) and 14 (western) samples from the two ordered lists new n 
were calculated for each sample by increasing n to 25 for each stratum. The final set of 50 
samples is listed along with their probabilities of inclusion in Table 1, and a map illustrating 
the spatial distribution of the samples is shown in Figure 1. In addition, a comparison of 
the forest type map used for sample selection with the forest types represented in our 
samples demonstrates that the design was effective at capturing the diversity of US forest 
types in the sample frame (Ruefenacht et al, 2008) (Figure 2). 
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2.2. Deriving the Scene-level Disturbance Products 


The creation of scene-level NAFD disturbance products has been described in detail 
previously (Huang et al, 2009; 2010a, 2010b, Thomas et al., 2011]. Here we briefly 
summarize the steps required to map disturbance within each Landsat sample scene and 
then describe the quality of the products in terms of severity and types of disturbance 
detected. 

At each of the 50 sample locations, a Landsat image time-series was constructed 
consisting of one image during the growing season (leaf-on conditions] for (initially] a 
target of every other year, between 1985 and 2005. Although biennial image acquisition 
was the initial targeted frequency, we were able to augment most image stacks with higher 
temporal frequency (i.e. mean intervals less than two years]. To optimize detection of 
change, images were chosen based on cloud cover and seasonality, with no attempt to 
synchronize skipped years in different stacks. Leaf-on seasonality required acquisition 
dates between June and September for most of the United States, although the range was 
extended to include May and October in the southern states. 

Each image stack was processed to maintain the highest radiometric and geometric 
standards (Huang et al., 2009]. Imagery was obtained as standard LIT (orthorectified at- 
sensor radiance] files from USGS EROS, and the latest version of the appropriate sensor 
calibration parameter set was applied. Geometric registration was checked and (as 
necessary] corrected by automated selection of image tie points and orthorectification 
(Huang et al., 2009]. The images were then converted to surface reflectance using the 
LEDAPS atmospheric correction package (Vermote et al., 1997; Masek et al., 2006] and 
assembled into a time series stack clipped to a common geographic extent. 
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Finally, water, cloud, and cloud shadow were identified and masked in each image. 
Water was mapped through a combination of decreasing reflectance with wavelength and 
low NDVI value (Huang et al., 2010a). Clouds were mapped using a set of visible/top-of- 
atmosphere temperature relationships (Huang et al., 2010b). Residual cloud contamination 
not identified in this step was also isolated as single-year "outliers" in the VCT forest 
disturbance analysis discussed below, and removed. 

Forest disturbances were mapped from the Landsat time series stacks using the VCT 
algorithm (Huang et al., 2010a). The algorithm used an automated approach to select forest 
training samples in each Landsat image and then calculated the distance in spectral space 
between each image pixel and the centroid of the forest training population (Huang et al., 
2008). Pixels close to the centroid of the forest population in the spectral space across the 
entire time period of record were classified as persistent forest for the entire observing 
period. Forest disturbance year for a given pixel was identified when that pixel's spectral 
properties exceeded an expected range of spectral deviation scores for at least two 
subsequent sequential time steps. The algorithm generated maps with classes for 
persistent non-forest, persistent forest, water, and the year of disturbance (Figure 3). 

These maps were filtered in both the spatial and temporal domains. Filtering was 
necessary to reduce false positive detections ("speckle") caused by residual image 
misregistration (Townshend et al., 2000; Knight & Lunetta, 2003). Disturbance pixel groups 
that were adjacent in space and time (i.e. within one time-period of each other) were 
identified. A majority filter was applied to achieve a minimum mapping (MMU) of 0.16 ha 
(two pixels) for static persistent forest and non-forest classes and 0.32 ha (4 pixels) for 
disturbance classes. Disturbance pixels in groups smaller than the MMU were converted to 
a majority class using decision rules in conjunction with a mode filter in the local 3x3 pixel 
neighborhood. Majority filtering typically reduced per-scene disturbed area by about 20%. 
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2.3. Per-scene Disturbance Product Quality 


Two approaches were used to describe disturbance product quality. First, six scenes 
representative of a range of forest and disturbance classes were selected to evaluate VCT 
performance (Thomas et al., 2011). These locations included a variety of disturbance types, 
including forest fire, harvest, thinning, land use conversion, and both storm and insect 
damage. The maps derived from VCT were compared to estimates determined 
independently using expert visual interpretation of Landsat and high-resolution satellite 
imagery. Because no field-based datasets exist that match the temporal density and spatial 
detail of the Landsat observation, a field-based validation was not undertaken (Cohen et ah, 
2010). Rather, trained interpreters, using established photo-interpretation guidelines, 
visually evaluated the Landsat TM imagery in tandem with high-resolution digital imagery 
such as that available on the Google Digital Earth. The Google imagery was used to aid 
determination of final land cover/use and to provide geographic context, while the 
determination of forest change came from visual interpretation of the Landsat TM imagery. 
The interpreters labeled change and no-change conditions using knowledge of both the land 
spectral properties and the spatial context of the landscape. This process was conducted at 
over random 600 points for each of the six validation scenes, with the points stratified by 
disturbance class (ie. year of disturbance mapped). For each scene, accuracy metrics of the 
VCT map were calculated by comparing that map with the visual interpretation results. 

The results yielded overall per-scene accuracies of 77% to 86% and kappa values of 
0.67 to 0.76 (Thomas et al., 2011). Considering only the forest change classes, omission and 
commission errors varied widely due to small sample sizes for individual year classes. 
Forest change users’ accuracies (100%-commission error%) averaged to 55-79% among 
the six scenes, although these accuracies increased by ~9% if the VCT year of disturbance 
was allowed to be within 1 year of the reference assignment. In general omission errors 
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were higher than commission errors. As a result, disturbance rates calculated from VCT 
map products were biased, underestimating total disturbance by an average of 24% across 
all validation sites. Omission errors mainly reflected the limitations of biennial sampling for 
capturing subtle disturbances (e.g. mechanical thinning], effects of seasonality within the 
image time series, and limitations of the VCT algorithm for mapping gradual, multi-year 
declines in forest health (i.e., stress] and change in the sparse forests of the western US. 

The specific type of disturbance was a lesser control on accuracy compared to the 
severity of disturbance, defined as the fraction of tree cover killed within the 0.32 ha MMU. 
Both natural and anthropogenic stand clearing disturbances (principally clear-cut harvest, 
severe fires, and major storm events] could be detected with high accuracy (75-85% 
detection accuracy], while non-stand clearing disturbances (including thinning, understory 
burns, and insect defoliation events] were only detected with 38% accuracy (Thomas et al., 
2011]. It should be noted that this latter figure increased to 60% accuracy if the allowable 
temporal window was relaxed to +/- 1 year. 


In a second, subsequent analysis (unpublished data] we used the TimeSync Landsat 
times series validation tool (Cohen et al. 2010] to better describe omission errors in terms 
of disturbance severity and type. TimeSync is a software environment to support 
visualization and interpretation of Landsat time series data. It includes simultaneous 
display of multitemporal Landsat image subsets, per-pixel Landsat spectral time series, and 
temporal snapshots of high-resolution imagery available in Google Earth. Interpreters can 


use the TimeSync environment to assign land cover, land use, and change labels to per-pixel 


trajectories, and then store this information in an online database. For this analysis we 
examined 150 plots randomly located within each of nine of the 50 sample scenes. Three 
severity classes were defined based on the percent of canopy removed or killed by 
disturbance (<34%, 34-67%, and >67%]. In the high and medium severity classes, 10% and 
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23% omission rates were observed. In contrast, the omission rate for the low severity 
disturbance class was 68%. However, while omissions error was linked to disturbance 
intensity, they were not uniform across disturbance type, indicating an interaction between 
severity and causal agent in terms of omission rates. 


For the three dominant types, harvest, fire and stress, omission rates in the low severity 
class were 76%, 38%, and 85%, respectively. This suggests that fire is the most likely and 
stress the least likely type of disturbance detected by VCT, with clear-cut harvest much 


more likely to be detected than selective harvest. For this analysis, we defined stress as an] 


observation of a spectral trend in the direction of disturbance that involved the loss of leaf 
area (or death] of live woody vegetation and was supported by observations within Google 


Earth. To be declared stress, the spectral trend had to be multi-year (typically 5 or more 


years]. We assumed that the large majority of such observations were associated with 


insect and disease activity (largely in conifer forests], recognizing that some observations 
were associated more directly with prolonged drought (e.g., in pinyon-juniper and 


ihrublands]. We could not observe seasonal losses of foliage associated with insects that 


would subsequently recover to full or near-full leaf area the following year (e.g., in eastern 
deciduous forests]. The TimeSync analysis supported the initial observation that omission 
error was greater than commission error and that omission error was more strongly 
associated with disturbance severity than type. 


2.4 Estimation of Area Disturbance Rates and Uncertainty 

To estimate disturbance rates, we first derived annual totals of disturbed area for each 
of the 50 NAFD sample scenes by interpolating across missing years of data within the 
scene. This was accomplished by allocating mapped disturbance from any multi-year 
period evenly to each year within that period. For example, if the 1992 image were missing 
from the time series, and 10,000 ha of disturbance were mapped during the 1991-1993 
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period, 5000 ha of disturbance would be allocated each to 1992 and 1993. Gaps in the 
image time series were asynchronous among samples, and more than half of the samples in 
each stratum provided data for any given year. Total forest area was also calculated for 
each scene by summing across three mapped classes: persistent forest, disturbed forest, and 
forest recovering from a pre-time-series disturbance (Huang et al., 2010a]. 

From our samples, total disturbed area per year and total forest area were estimated for 
the two strata using the generalized Horvitz-Thompson estimator for unequal-probability 
designs (Horvitz & Thompson, 1952]: 


■ 4=2 


A 


[ 1 ] 


where A is the unbiased estimate of either forest area or disturbed area for a given year, y, is 
the forest or disturbance area for the i th sample in that year, and m is the probability of 
inclusion for the i th sample. Proportional forest disturbance rate by stratum was calculated 
by dividing this total by the estimated total forest area in the stratum. Estimation of 
variance was achieved using the conservative collapsed stratum variance method (Cochran, 
1977], 


3.0 Results 


3.1. Per-Scene Disturbance Rates 

Raw, per-scene forest area and disturbance rates for each of the 50 samples showed 
significant geographic variability (Figure 4], Disturbance rates in the southern United 
States and Pacific Northwest were consistently high, both in terms of absolute area and 
fraction of mapped forest cover affected. Samples from the Carolinas, the Gulf Coast, and 
Northern Louisiana exhibited fractional disturbance rates in the range of 1.5 - 2.0%/yr. 
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Absolute and fractional disturbance rates in the northeast were considerably lower (< 
0.5%/yr]. While the absolute area of disturbance in the Interior (mountain] West was low, 
the forest cover base in these locations is limited and generally confined to mountain 
environments, and corresponding fractional disturbance rates were quite variable. 

3. 2. Eastern versus Western US Disturbance 

Having selected independent sample sets for the eastern and western United States, 
disturbance rates and variances were estimated separately for each stratum (Figure 5]. We 
found similar average disturbance rates (1.10%/yr and 1.05%/yr, respectively] for eastern 
and western forests, and significant interannual variability in disturbance rates within each 
of the stratum. However, the annual rates of disturbance were asynchronous in the two 
strata. Although both strata showed decreases in disturbance rate from the late 1980's to 
the early 1990's, thereafter the trajectories within the two regions diverged, with 
disturbance rates peaking in the east in the late 1990s and in the west in the early 2000s. 

The national sampling approach does not support disaggregating the results to finer 
geographic areas since the probability of inclusion for each sample scene was calculated 
assuming two strata. Thus there is no way to produce a subregional area estimate other 
than simple averaging (ie. assumed equal probability of inclusion]. With that caveat, we can 
calculate such "apparent" disturbance rates for smaller groups of scenes to interrogate 
regional patterns. Disturbance rates were explored independently for groups of scenes 
within the northern and southern subregions of the eastern stratum, and for groups of 
scenes within the intermountain and coastal subregions of the western stratum (Figure 6]. 
These boundaries correspond to the US Forest Service FIA regions. In the western US, 
disturbance rates in the intermountain west were generally lower than the coastal region, 
with the difference becoming more pronounced in the late-1990s. In the eastern US, the 
northern and southern subregions show markedly different temporal patterns. While the 
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northern subregion indicates a stable, low rate of disturbance, disturbance rates in the 
southeast are much higher and vary significantly on 5-10 year timescales. It should be 
noted, however, that because the results for individual subregions do not come from the full 
per-stratum sample, variances cannot be calculated for these rates. Thus while the 
temporal patterns can be broadly interpreted, their statistical significance is uncertain. 

3.3 National Disturbance Rates 

Aggregating the results from the eastern and western strata gives an estimate of 
national forest disturbance rates since 1985 (Figure 7). We find that an average of 2.77 
Mha/yr +/- 0.36Mha [la] of US forestland was disturbed each year during the 1985-2005 
period. This figure corresponds to a fractional disturbance rate of 1.09% forestland per 
year. The net underestimation bias (see Section 2.3) suggests that the VCT algorithm 
missed about 24% of actual disturbance that can be detected visually within the Landsat 
imagery. Adjusting the NAFD results for this underestimation bias would increase the 
national disturbance rate from 2.77 Mha/yr (1.09% forestland/yr) to 3.63 Mha/yr (1.42% 
forestland/yr). 

Unlike studies that report only average disturbance rate over multi-year periods our 
results suggest that forest disturbance rate can change rapidly across large regions (Figure 
7). The time series of national disturbance varies about the mean value, from a low 0.8% of 
US forest area in 1992 to a high of 1.4% in 2000. In general, values during the early 1990s 
were consistently lower than average, while values in the late 1990s were consistently 
higher. The uncertainty bounds based only on sampling error are less than the magnitude 
of the interannual variability, suggesting that these swings in disturbance rates at the 
national scale are "real" and not simply artifacts of the sampling. However, some caution 
must be taken given that only sampling errors are shown on Figures 5 and 7; inclusion of 
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measurement errors would increase estimated error and reduce confidence in the 
estimates of interannual variability. 

3.4. US Forest Area 

Total forest area for the conterminous US was also calculated from the NAFD results, 
estimated as the sum of permanent (undisturbed] forest and all annual disturbance classes 
within the VCT classification. The estimated total forest areas for the eastern and western 
strata were 182 x 10 6 ha and 74 x 10 6 ha, respectively, for a total estimate of 255 x 10 6 ha. 
The US Forest Service reported areas of 157 x 10 6 ha and 96 x 10 6 ha for eastern and 
western forestland in 2007, respectively (Smith et al., 2009]. Thus while the NAFD national 
forest area corresponds closely to the US Forest Service national estimate of 253 x 10 6 ha 
for forestland in the conterminous US, the per-stratum areas are somewhat different. 


4. Discussion 

4.1. National Disturbance Rates 

The results presented here represent the first satellite-based estimates of yearly forest 
disturbance for the conterminous United States. From our geographic sample, we have 
estimated that, on average, 1.09% of US forestland was affected by disturbances each year 
during the 1985-2005 epoch. It should be noted, however, that our sample was primarily 
based on representing forest type, and the geographic distribution of disturbances does not 
necessarily correspond to that of forest types. Therefore there is some additional 
uncertainty in the national rates not captured by the estimate of sampling error. 

It is of interest to compare our results to data from US land management inventories, as 
well as previous remote sensing studies (Table 2]. As noted in the Introduction, there is no 
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single database from the US forest inventory that tracks the area of forest disturbance. 
However, there are independent reports for the area of harvest (including area affected by 
partial harvest) (Smith et al., 2009), the area affected by forest fire (EPA, 2011), and the 
forest area affected by insect mortality (USDA Forest Service, 2010). The age distribution of 
US forests reported in the FIA can also be converted to an average stand-clearing 
disturbance rate assuming that all forestland area younger than 20-years experienced a 
stand-clearing disturbance event during the last 20 years. This calculation implicitly 
excludes non stand-clearing events that would reduce stand biomass but not necessarily 
alter the FIA age assignment. 

Previous remote sensing estimates of US forest disturbance are given in Masek et al 
(2008) based on a 10-year (1990-2000) wall-to-wall analysis of Landsat imagery, and 
Hansen et al (2010) based on 500-m resolution MODIS imagery for the 2000-2005 time 
period. Like the NAFD sample, these estimates represent gross forest cover loss; that is 
they include both losses due to disturbance (but that will return to forest cover) as well as 
losses due to permanent forest conversion (deforestation). However, the rate of 
deforestation in the US is relatively low compared to the turnover due to harvest and other 
forms of disturbance (Smith et al., 2009). Drummond and Loveland (2010) also used 
Landsat data in 5-7-year epochs to estimate land cover transitions in the eastern US. They 
found gross forest conversion (which includes harvest and permanent conversion) of 0.34 
Mha/yr, although this cannot be directly related to the NAFD eastern stratum since the 
Drummond and Loveland (2010) study area only included areas east of the Mississippi and 
Ohio Rivers. 

Comparing these studies suggests a general consistency among the satellite-based 
disturbance estimates (Table 2). The NAFD national results are higher than either Masek et 
al (2008) or Hansen et al (2010), presumably reflecting an improved ability to identify 
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subtle disturbance signals when using finer spatial resolution data (compared to MODIS] 
and shorter mapping intervals (compared to the decadal Landsat study). The NAFD rates 
are also higher than those calculated from the FIA forestland age distribution, suggesting 
that the NAFD results are identifying a significant amount of partial harvest (i.e., thinning) 
that would not be severe enough to reset the FIA stand age. This is corroborated by the 
TimeSync assessment that indicated most moderate and some low severity disturbances 
were detected. 

All of the satellite-based estimates are considerably lower than the rates obtained 
by summing the inventory-based estimates of forest harvest, fire, and insect mortality. The 
total area of disturbance obtained from the inventory reports is about 8.1 Mha/yr, with 4.3 
Mha/yr attributable to harvest activity alone (Table 2). There are two likely explanations 
for the discrepancy. First, the inventory data generally record the area affected by a type of 
disturbance rather than the actual area of forest disturbed. For example, the insect 
mortality estimate is based on the area enclosed by areal sketch map polygons, while only a 
fraction of the stands within those polygons have succumbed (Johnson and Wittwer 2008). 
Meddens et al (2012) assessed the actual tree mortality area within USFS Aerial Detection 
Survey polygons for beetle outbreaks in the western US, and estimated up to 380,000 ha/yr 
of actual tree mortality averaged for the 1997-2010 period. Similarly the area considered 
affected by partial harvest (e.g. selective cutting or strip cuts) may be less than the actual 
area cleared. In contrast the satellite-based results offer the potential for a more exact 
estimate of disturbed crown area. 

However, some of the discrepancy also likely reflects limitations of the current 
NAFD methodology. The methods described here are sensitive to most severe and moderate 
disturbances, but commonly miss events that kill a small fraction of trees within the stand 
or disturbances that emerge gradually (such as drought or insect mortality). In addition, 
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the accuracy of the VCT change detection was hampered by occasional inconsistencies in 
the acquisition date of images within the time series. Using the bias-adjusted figure derived 
from Thomas et al. (2011) of 3.63 Mha/yr (1.42% forestland/yr) brings the satellite-based 
disturbance estimates closer to the inventory figure. Finally, there remains the possibility 
that the NAFD sampling approach introduced some bias into the results compared to 
geographically comprehensive estimates from forest inventories. 

4.2. Interannual Variability 

The NAFD results also suggest significant interannual changes in disturbance rates. 
However, our national estimates integrate across the many factors (forest policy and 
management, fire ignition, local economic conditions, land ownership) that determine forest 
disturbance rates when viewed at the local or regional scales. We do know that a primary 
causal factor in disturbance rates is forest management, specifically harvest. The quantity 
of timber harvested is a function of societal, economic and political processes (Prestemon 
2001). Harvest rates are heavily influenced by volatile wood product markets which feed 
into and are fed by overall economic activity as measured by GDP (Daniels et al, 
2005;Luppold and Miller, 2005 

In the western stratum the NAFD results show an early decrease and then significant 
increases, peaking in 2000. These dynamics likely relate to decreased western harvests in 
the early 1990's, increased forest fires, particularly in the interior west in the late 1990's 
(Daniels, 2005; mtbs.gov), and increased insect activity (Meddens et al., 2012) (Figure 8a). 
The changes in fire dynamics follow the widely recognized pattern of drier and warmer 
conditions, more frequent fires, and more intense insect outbreaks over the last two 
decades in the western US and Canada (Williams et al. 2010; Westerling et al., 2006; Kurz et 
al, 2008). 
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In the eastern US, forest harvest appears to the dominant mechanism for disturbance, 
particularly is the southeastern US suggesting that these rates may be related to US 
economic activity as reflected in the Gross Domestic Product (GPD] which strongly reflects 
changes in housing starts (Figure 8b] (Napton et al, 2010]. Recessions of 1990 to 1991 and 
2000 to 2001 caused drops in the rate of increase of GDP that were mirrored by drops in 
eastern forest disturbance, and the period of highest disturbance rate (1997 to 2000] was 
also a period of large year over year increase in GDP. 

Some caution must be exercised in interpreting these trends however. A detailed 
analysis of the per-scene results indicates that much of the interannual variability in the 
eastern stratum was driven by accelerated harvest between 1996 and 2000 in two sample 
scenes in Oklahoma (WRS-2 p26r36] and eastern Texas (WRS-2 p26r37] (Schleeweis, 
2012]. The sudden increase in harvest rates was noted at the time by the Assistant State 
Forester for Oklahoma (K. Atkinson, pers. comm], and was attributed to a combination of 
regional timber prices, local increases in demand and capacity due to the opening of new 
chip and stud mills, and real estate exchanges by timber companies. 

While the scene-level Landsat VCT products successfully captured this local acceleration 
in disturbance, such isolated (in space and time] jumps in disturbance rate were not built 
into the original sampling framework. The two of the four drivers of the sampling scheme 
were the capture of a range of forest types (forest diversity] and the maximization of 
sample geographic dispersion. Essentially, the existence of isolated "pockets” of rapid 
disturbance, which could not have been predicted a priori and thus could not be 
incorporated into the sampling framework, raises questions as to the efficacy of using a 
limited geographic sample to fully capture variability in disturbance rates. Similar points 
have been made regarding tropical deforestation, which also exhibits a geographic 
"clumping" that is difficult to predict a priori (Tucker and Townshend, 2000]. 
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4.3. Sample size and sampling error 


The analysis of the 50-scene NAFD sample was conducted in two phases: a first phase 
using 23 scenes, and a subsequent phase with an additional 27 scenes. This gives us an 
opportunity to examine the effect of increased sample size on stratum-level estimates. Mean 
disturbance rate and variances were calculated separately for the samples from the first 
phase and the entire 50-scene set. Adding more scenes to the sample should improve the 
precision of the disturbance estimates, and the new rates should typically he within two 
standard deviations of the estimates from the smaller set. 

Comparing the Phase 1 and combined Phase 1 and 2 results (Figure 9] for the eastern 
stratum, the overall temporal pattern is similar, and the individual rates from the full set 
actually lie within the la of the smaller set. However, in the western stratum, estimates of 
rates during particular eras (mid-1980's, 2000-2002) are separated by at least two 
standard deviations. While it is reasonable to assume that the full set of 50 scenes provides 
a more precise estimate than the smaller Phase 1 sample, this again illustrates the difficulty 
in using a limited geographic sample to capture regional disturbance rates. 

5.0. Conclusions 

The NAFD project has used a geographic sample of Landsat scenes to quantify forest 
disturbance rates across the conterminous United States for the 1985-2005 epoch. We 
found that an average of 2.77 Mha/yr (1.09% forestland/yr) was disturbed each year 
during this epoch. Adjusting this figure for the net measurement error (bias) would 
increase it to 3.63 Mha/yr (1.42% forestland/yr). Both of these figures are lower than the 
combined area affected by harvest, fire, and insect damage reported from land management 
inventories, although the figure is higher than the disturbance rate of forestland expected 
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from the FIA age distribution. Much of the discrepancy is associated with low severity 
disturbances in which a small fraction of the canopy cover is killed. 

Disturbance rates in the western US appear to have increased substantially 
beginning in the late 1990s and closely track increased fire frequency in the intermountain 
west during this period. Rates in the eastern US more closely tracked variability in 
economic activity (GDP], presumably reflecting fluctuations in harvest. The temporal 
variability in eastern disturbance was greater than expected, although it is unclear whether 
these variations reflect national patterns or local "pockets" of increased disturbance. 

The work presented in this study points the way toward a consistent monitoring 
framework for US forest dynamics based on interannual satellite observations. Disturbance 
rates need to be quantified at the spatial resolution where human management and natural 
disturbances occur and at a time step relevant to driving processes. Although available 
forest inventories provide invaluable information on US forest dynamics, they have not 
been designed to provide a consistent geospatial view of disturbance. As a result, it is not 
clear from those data how long-term changes in land management and natural disturbance 
rates ultimately affect ecosystem structure and functioning across the nation. In principle, 
satellite observations can help meet this need. 

The NAFD results described here relied on a limited geographic sampling, and the 
divergence in results between the initial set of 23 samples and the full set of 50 highlights 
the difficulty in using such a sampling approach to constrain disturbance dynamics. Our 
detailed analysis of the regional and local variations in the observed disturbance rates has 
revealed that the geography of processes that drive disturbance and specific local 
disturbance events have a strong influence over variations in area-averaged disturbance 
rates from year to year (Schleweiss. 2012]. As a result we have turned toward wall-to-wall, 
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annual mapping in the current phase of the NAFD project to circumvent these limitations. 
The advent of free Landsat data distribution, coupled with increased computing power, has 
led to the use of "every clear pixel" in mapping land dynamics (Wulder et al., 2012; Zhu et 
al., 2012]. Availability of such comprehensive information will precipitate more advanced 
ecosystem process models that can ingest the richness of these data, and lead to a much 
improved ability to monitor and forecast ecosystem responses to climate change and 
human management. Coupled with repeat-measurement inventory data, when those 
become available over the next several years, we will likely soon have a powerful system for 
spatially explicit monitoring of forest cover in the conterminous US. 
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Table Captions 


Table 1. Ordered list of the sample scenes chosen in each stratum (eastern, western) with 
probability of inclusion. 


Table 2. Comparisons of estimated US disturbance rate from this study with previous 
satellite and inventory-based studies 
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Table 1. Ordered list of the sample scenes chosen in each stratum 
(eastern, western) with probability of inclusion. 


Eastern stratum Western stratum 


WRS-2 Probability of WRS-2 Probability of 


Order 

path/row 

inclusion (n) 

path/row 

inclusion (n) 

1 

21/37 

0.138 

37/34 

0.117 

2 

27/38 

0.130 

35/34 

0.234 

3 

22/28 

0.154 

34/37 

0.183 

4 

18/35 

0.146 

37/32 

0.157 

5 

25/29 

0.118 

47/28 

0.183 

6 

17/31 

0.142 

36/37 

0.213 

7 

19/39 

0.126 

41/32 

0.213 

8 

16/36 

0.177 

43/33 

0.198 

9 

26/36 

0.157 

41/29 

0.162 

10 

12/31 

0.122 

35/32 

0.112 

11 

21/39 

0.134 

45/29 

0.173 

12 

nm 

0.146 

44/26 

0.137 

13 

26/34 

0.142 

42/35 

0.137 

14 

24/37 

0.114 

42/28 

0.208 

15 

21/30 

0.197 

44/29 

0.168 

16 

23/35 

0.130 

46/32 

0.188 

17 

26/37 

0.114 

46/31 

0.223 

18 

14/31 

0.142 

42/29 

0.157 

19 

16/35 

0.142 

46/30 

0.178 

20 

16/41 

0.169 

47/27 

0.157 

21 

23/28 

0.165 

48/27 

0.168 

22 

20/33 

0.154 

40/37 

0.152 

23 

15/31 

0.122 

34/34 

0.193 

24 

27/27 

0.087 

33/30 

0.183 

25 

19/36 

0.122 

45/30 

0.152 
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Table 2. Comparisons of estimated US disturbance rate from this study with previous 
satellite and inventory-based studies 


Category 

Source 

Period of 
record 

Extent 

Method/source 

Type 

Forest area 

disturbed 

(Mha/yr) 

Remote 

Sensing 

Masek et 
al. (2008) 

1990-2000 

Conterminous 

US 

Wall-to-wall 
mapping with 
two dates of 
Landsat 
imagery 

Observed, 

stand- 

replacing 

disturbances 

1.97 


Hansen et 
al. (2010) 

2000-2005 

US (incl. AK 
and HI) 

Sample-based 
mapping using 
two dates of 
Landsat 
imagery and 
MODIS 

Observed 

disturbance 

2.5 


This study 

1985-2005 

Conterminous 

US 

Sample-based 
mapping with 
dense time 
series 

Observed 

disturbance 

2.77 (+/- 
0.76) 






Adjusted 

disturbance 1 

3.63 

FIA Age 


1988-2008 

Conterminous 

US 

Annualized 
area of forests 
<20 years in 
age 

Stand- 

replacing 

disturbance 

2.2 

Inventory 

Smith et al. 
(2009) 

2001-2005 

Conterminous 

US 

Remeasured 
inventory plots 
in east, 

National Forest 
activity reports 
in west 

Clearcut and 

Partial 

harvest 

4.34 


US EPA 
(2010) 

1990-2005 

Conterminous 

US 


Forest fire 

1.01 


USD A 
Forest 
Service 
(2010) 

Late- 
2000 ’s 

US (incl. AK, 
HI) 

Aerial 
Detection 
Survey (ADS) 

Area 

affected by 

insect 

disturbance 

2.74 


Meddens et 
al. (2012) 

1997-2010 

Western US 

ADS calibrated 
to area of 
mortality 

Beetle 

mortality 

0.38 


1 Adjusted disturbance includes lower-intensity forest disturbance and is based on application of the 
bias omission factor as described in Section 4.1. 
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Figure Captions 


Fig. 1. NAFD sample scenes (unique, non-overlapping scene areas) overlaid on the US 
Forest type map (Ruefenacht et al., 2008). 


Fig. 2: Fraction of forest types captured by the NAFD sample for the (a) eastern and (b) 
western strata (grey bars) compared to the actual area of each type from the map of 
Ruefenacht et al. (2008) (black bars). 

Fig. 3. Example of Vegetation Change Tracker (VCT) algorithm output for Landsat WRS-2 
path 43 row 33 (Sierra Nevada, California). The VCT uses per-pixel annual Landsat time 
series to generate maps with classes for permanent forest, non-forest, and disturbance 
occuring in each year. Disturbed patches are color-coded by year of disturbance. Lake 
Tahoe is depicted as the blue area in the north-center of the full scene. The zoomed image 
shows small patches of harvest as well as the large burn scar of the 1992 Cleveland fire. 

Lig. 4: Per-scene results of the NAFD disturbance analysis. The overall length of the bar 
shows total NAFD derived forest land area for the WRS-2 sample location (x-axis). The black 
portion of the bar is equal to the area of forest mapped as disturbed during 1985-2005. 

Note that samples along coastlines show low forest area due to the large proportion of 
water that falls within the samples footprint. 

Fig. 5: Disturbance rates and uncertainty estimates calculated separately for the eastern 
(dark grey) and western (light grey) strata. Error bars represent sampling error (+/- la) 
and western points are offset by 0.2 years for clarity. 
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Fig. 6. Disturbance rates (% forest change per year) calculated for (a) the south and north 
subregions of the eastern stratum and (b) the intermountain and coastal subregions of the 
western stratum. The thickness of each shaded area gives the disturbance for that 
subregion, such that the upper surface gives the total disturbed area for each stratum. 

Fig. 7. NAFD national disturbance rates obtained by combining results for the eastern and 
western strata, calculated by year (% of conterminuous forest area disturbed per year). 
Error bars represent sampling error (+/- la). 


Figure 8. Temporal patterns of forest disturbance rates compared to factors 
expected to influence those rates for the western (a) and eastern (b) U.S. Wood 
production and GDP rates refer to year-to-year changes in 2-year running average 
national estimates (Howard, 2007). Fires are areal estimates of forest fire 
disturbance alone for the western U.S.obtained from the Monitoring Trends in 
Burn Severity (MTBS) project (mtbs.gov). Insects represent the upper bound 
of beetle-associated mortality area (Meddens et al., 2012). 

Fig. 9. Comparison of Phase I (23 scenes, top) and Phase II (50 scenes, bottom) estimated 
disturbance rates (percent forest cover per year) and sampling error for the eastern (left) 
and western (right) strata. 
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Other western softwoods 
Pinyon/Juniper group 
Ponderosa pine group 
Redwood group 
Tanoak/iaurel group 
Western larch group 
Western oak group 
Western white pine group 


Aspen/bi rch group 
| Elm/ash/cottonwood group 
| Loblolly/Shortleaf pine group 
| Longleaf/Slash pine group 
Oak/Hickory group 
Oak/Pine group 


Oak/gum/Cypress group 
| Spruce/Fir group 
| Tropical hardwoods 
| White/red/jack pine group 
| Maple/beech/birch group 


□ NAFD samples 

FIA forest group type 

Non-forest 
Alder/maple group 
1 California mixed conifer 
| Douglas-fir group 
| Fir/spruce/mountain hemlock 
| Hemlock/sitka spruce group 
| Lodgepole pine group 
■ Other western hardwoods 


Fig. 1. NAFD sample scenes (unique , non-overlapping scene areas) overlaid on the US 
Forest type map (Ruefenacht et al., 2008). 
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Fig. 2: Fraction of forest types captured by the NAFD sample for the (a) eastern and (b] 
western strata (grey bars) compared to the actual area of each type from the map of 
Ruefenacht et al. (2008) (black bars). 


39 





Non-forest 
Perm. Forest 
Water 
<1984 

1986 

1988 

1990 

1992 

1994 

1996 

1998 

2000 

2002 

2004 


4 km 


Fig. 3. Example of Vegetation Change Tracker (VCT] algorithm output for Landsat WRS-2 
path 43 row 33 (Sierra Nevada, California). The VCT uses per-pixel annual Landsat time 
series to generate maps with classes for permanent forest, non-forest, and disturbance 
occuring in each year. Disturbed patches are color-coded by year of disturbance. Lake 
Tahoe is depicted as the blue area in the north-center of the full scene. The zoomed image 
shows small patches of harvest as well as the large burn scar of the 1992 Cleveland fire. 
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Fig. 4: Per-scene results of the NAFD disturbance analysis. The overall length of the bar 
shows total NAFD derived forest land area for the WRS-2 sample location (x-axis). The black 
portion of the bar is equal to the area of forest mapped as disturbed during 1985-2005. 

Note that samples along coastlines show low forest area due to the large proportion of 
water that falls within the samples footprint. 
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Year 


Fig. 5: Disturbance rates and uncertainty estimates calculated separately for the eastern 
(dark grey) and western (light grey) strata. Error bars represent sampling error (+/- la) 
and western points are offset by 0.2 years for clarity. 
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Fig. 6. Disturbance rates (% forest change per year] calculated for [a] the south and north 
subregions of the eastern stratum and (b] the intermountain and coastal subregions of the 
western stratum. The thickness of each shaded area gives the disturbance for that 
subregion, such that the upper surface gives the total disturbed area for each stratum. 
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Fig. 7. NAFD national disturbance rates obtained by combining results for the eastern and 
western strata, calculated by year (% of conterminuous forest area disturbed per year). 
Error bars represent sampling error (+/- la). 


44 



2.0 


a) Western forests 


b) Eastern forests 


40 



Figure 8. Temporal patterns of forest disturbance rates compared to factors 
expected to influence those rates for the western (a] and eastern (b] U.S. Wood 
production and GDP rates refer to year-to-year changes in 2-year running average 
national estimates (Howard, 2007]. Fires are areal estimates of forest fire 
disturbance alone for the western U.S.obtained from the Monitoring Trends in 
Burn Severity (MTBS] project (mtbs.gov]. Insects represent the upper bound 
of beetle-associated mortality area (Meddens et ah, 2012). 
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Fig. 9. Comparison of Phase I (23 scenes, top] and Phase II (50 scenes, bottom] estimated 
disturbance rates (percent forest cover per year] and sampling error for the eastern (left] 
and western (right] strata. 
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